function rtp = elements2rtp(el, ta, mu)

r = el.a*(1-el.e^2) / (1+el.e*cos(ta));
%[km] position magnitude from classical orbital elements to rtp coordinate
%frame

rrth = [r; 0; 0];
%[km] position vector in rth coordinate frame

v = sqrt(mu*(2/r - 1/el.a));
%[] speed of the space craft

gamma = atan(r*el.e*sin(ta) / (el.a*(1-el.e^2)));
%[] Flight path angle

vrth = [v*sin(gamma); v*cos(gamma); 0];
%[] velocity vector in rth coordinate system

theta = ta+el.aop;
%[] Theta for transformation matrix

Arth2xyz = [cos(el.raan)*cos(theta)-sin(el.raan)*cos(el.i)*sin(theta), -cos(el.raan)*sin(theta)-sin(el.raan)*cos(el.i)*cos(theta), sin(el.raan)*sin(el.i);
            sin(el.raan)*cos(theta)+cos(el.raan)*cos(el.i)*sin(theta), -sin(el.raan)*sin(theta)+cos(el.raan)*cos(el.i)*cos(theta), -cos(el.raan)*sin(el.i);
            sin(el.i)*sin(theta), sin(el.i)*cos(theta), cos(el.i)];
%[] Transformation matric from rth to xyz

rxyz = Arth2xyz * rrth;
%[] rotate position vector from rth to xyz coordinate frame

vxyz = Arth2xyz * vrth;
%[] rotate velocity vector from rth to xyz

theta = acos(rxyz(3)/r);
%[] angle theta measured from z axis down to position vector.

phi = atan2(rxyz(2), rxyz(1));
%[] angle phi measured from positive x axis to the projection of the
%potision vector r onto the xy plane

rtp = [r; theta; phi];
%[] Position of the spacecraft with r theta phi.

Axyz2rtp = [sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta);
    cos(theta)*cos(phi), cos(theta)*sin(phi), -sin(theta);
    -sin(phi), cos(phi), 0];
%[] rotation matrix from xyz to rtp

vrtp = Axyz2rtp * vxyz;
%[] Velocity of the spacecraft in rtp frame.

rtp = [rtp' vrtp']';
%[] state of the spacecraft in rtp coordinate frame.

end